#ifndef BTL_FEATURES_FEATUREDESCRIPTOR
#define BTL_FEATURES_FEATUREDESCRIPTOR

#include <btl/Maths/OnlineStatistics.hpp>
#include <Eigen/Dense>

namespace btl
{
namespace features
{

class FeatureDescriptorBase
{

   public:

      virtual ~FeatureDescriptorBase();

};

template <class VectorValueType>
class FeatureDescriptor : public FeatureDescriptorBase
{

   public:

      typedef VectorValueType value_type;
      typedef Eigen::Map< Eigen::Matrix<value_type, Eigen::Dynamic, 1> > descriptor_vector;

      virtual ~FeatureDescriptor(void);

      FeatureDescriptor& operator=(const FeatureDescriptor&);

      virtual descriptor_vector asVector() = 0;

};

template <typename PixelType, int Width, int Height>
class PatchFeatureDescriptor : public FeatureDescriptor<PixelType>
{

   public:

      typedef typename FeatureDescriptor<PixelType>::descriptor_vector descriptor_vector;
      typedef Eigen::Matrix<double, PixelTraits<PixelType>::NUM_CHANNELS, 1> pixel_vector;

      PatchFeatureDescriptor(void);

      PatchFeatureDescriptor(const ImageRegion<const PixelType>& image, int x, int y);

      virtual ~PatchFeatureDescriptor(void);

      int width(void) const;

      int height(void) const;

      ImageRegion<const PixelType> patch(void) const;

      pixel_vector mean(void) const;
      pixel_vector standardDeviation(void) const;

      descriptor_vector asVector();

   protected:
      void calculateStats(void);

      SmallImage<PixelType, Width, Height> _patch;

      pixel_vector _mean;
      pixel_vector _stdDev;

};

template <int Width, int Height>
double ncc(
   const PatchFeatureDescriptor<PixelGrey, Width, Height>& fd1,
   const PatchFeatureDescriptor<PixelGrey, Width, Height>& fd2);

// other descriptors to implement
// (note: some will use OpenCV)
// - SURF
// - SIFT
// - Denis' descriptor
// - Haar wavelet coefficients

} // namespace features
} // namespace btl

namespace btl
{

using features::FeatureDescriptor;
using features::PatchFeatureDescriptor;

} // namespace btl


// ====================================================================
// === Implementation


namespace btl
{
namespace features
{

inline FeatureDescriptorBase::~FeatureDescriptorBase() {}

template <class VectorValueType>
inline FeatureDescriptor<VectorValueType>::~FeatureDescriptor(void) {}

template <class VectorValueType>
inline FeatureDescriptor<VectorValueType>& FeatureDescriptor<VectorValueType>::operator=(const FeatureDescriptor&)
{
   return *this;
}

template <typename PixelType, int Width, int Height>
inline PatchFeatureDescriptor<PixelType, Width, Height>::PatchFeatureDescriptor(void) {}

template <typename PixelType, int Width, int Height>
inline PatchFeatureDescriptor<PixelType, Width, Height>::PatchFeatureDescriptor(const ImageRegion<const PixelType>& image, int x, int y)
{
   const int halfW = Width / 2;
   const int halfH = Height / 2;
   _patch.copy(image.region(x - halfW, y - halfH, Width, Height));

   this->calculateStats();
}

template <typename PixelType, int Width, int Height>
inline PatchFeatureDescriptor<PixelType, Width, Height>::~PatchFeatureDescriptor(void) {}

template <typename PixelType, int Width, int Height>
inline int PatchFeatureDescriptor<PixelType, Width, Height>::width(void) const
{
   return _patch.width();
}

template <typename PixelType, int Width, int Height>
inline int PatchFeatureDescriptor<PixelType, Width, Height>::height(void) const
{
   return _patch.height();
}

template <typename PixelType, int Width, int Height>
inline ImageRegion<const PixelType> PatchFeatureDescriptor<PixelType, Width, Height>::patch(void) const
{
   return _patch;
}

template <typename PixelType, int Width, int Height>
inline void PatchFeatureDescriptor<PixelType, Width, Height>::calculateStats(void)
{
   typedef PixelTraits<PixelType> traits_type;
   const int NUM_CHANNELS = traits_type::NUM_CHANNELS;
   btl::maths::OnlineStatistics stats[NUM_CHANNELS];
   for(int i = 0; i < _patch.height(); ++i)
      for(int j = 0; j < _patch.width(); ++j)
         for(int k = 0; k < NUM_CHANNELS; ++k)
            stats[k].push(double(traits_type::get(k, _patch[i][j])));

   for(int k = 0; k < NUM_CHANNELS; ++k)
   {
      _mean(k) = stats[k].mean();
      _stdDev(k) = stats[k].populationStdDev();
   }
}

template <typename PixelType, int Width, int Height>
inline typename PatchFeatureDescriptor<PixelType, Width, Height>::pixel_vector
PatchFeatureDescriptor<PixelType, Width, Height>::mean() const
{
   return _mean;
}

template <typename PixelType, int Width, int Height>
inline typename PatchFeatureDescriptor<PixelType, Width, Height>::pixel_vector
PatchFeatureDescriptor<PixelType, Width, Height>::standardDeviation() const
{
   return _stdDev;
}

template <typename PixelType, int Width, int Height>
inline typename PatchFeatureDescriptor<PixelType, Width, Height>::descriptor_vector PatchFeatureDescriptor<PixelType, Width, Height>::asVector()
{
   return descriptor_vector(_patch.data(), _patch.width() * _patch.height());
}

template <typename FeatureDescriptorType>
inline double l2norm(const FeatureDescriptorType& fd1, const FeatureDescriptorType& fd2)
{
   // FIXME: implement
   return 0.0;
}

inline double ssd(const FeatureDescriptor<unsigned char>& fd1, const FeatureDescriptor<unsigned char>& fd2)
{
   return 0.0;
}

template <int Width, int Height>
inline double ncc(const PatchFeatureDescriptor<PixelGrey, Width, Height>& fd1, const PatchFeatureDescriptor<PixelGrey, Width, Height>& fd2)
{
   double ncc = 0.0;
   double normer = 1.0 / (fd1.standardDeviation()(0) * fd2.standardDeviation()(0));

   const double meanA = fd1.mean()(0);
   const double meanB = fd2.mean()(0);

   for(int i = 0; i < Height; ++i)
   {
      for(int j = 0; j < Width; ++j)
      {
         double ax = fd1.patch()[i][j] - meanA;
         double bx = fd2.patch()[i][j] - meanB;
         ncc += ax * bx * normer;
      }
   }

   ncc /= Width * Height;

   return ncc;
}


} // namespace features
} // namespace btl


#endif //BTL_FEATURES_FEATUREDESCRIPTOR

